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I. Introduction 

The equations governing coapreiaible vlacoua flow have been known for 
■ore than a century < Although aany apeclal-case solutions have been 
determined amlytically over the years, aany others of Intereat have 
continued to defy ■athenatlcal analysis. Wind tunnels and other experi- 
mental facilities have served as invaluable tools in the integration, by 
physical slasulatlon, of these equations of fluid Botion. During the 
last decade the computer has come to share — through its use of numeri- 
cal simulations - the work of the earlier analytical and experimental 
tools in determining new flow solutions. 

Like the limits on the range of problems that can be solved analytically, 
there are limits on the range of flow cases that can be accurately simu- 
lated in experimental facilities. The experimental limits are imposed 
by such factors as tunnel size, wall interference, and stream uniform- 
ity [1]. Similarly, the range of computer flow simulations is also 
limited, principally by computer speed and memory storage. Fortunately, 
the limits of the theoretical, experimental, and computational techniques 
are different; as a result, the range of applicability afforded by the 
three is greater than that attainable with any one of them. Moreover, 
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In rtglonc whtrt they overlap, one epproech cen be uecd to verify the 
reeulte of another. 

Ncvcrthelees , ve atlll cannot aolve ell fluid flow problOBa of intercat, 
nor can we anticipate that capability in the near future. However, 
becauee of the present rapid and potential large growth of couputer 
capabilities, Buch eophasls ia being placed on the developBent of cob- 
pucational fluid dynamics [2,3], Before we can calculate the flow field 
about a complete aircraft configuration at flight Reynolds numbers, 
there will have to be great progress in developing powerful and reliable 
computer hardware, in understanding md modeling the physics of turbulent 
flow, and in devising accurate and efficient nuaerical methods. That 
progress will depend, to a significant degree, on theoretical, experi- 
mental, and numerical research. f)ne clement, the devising of an effi- 
cient numerical method, is discussed in this paper. 

During the last 10 years many significant contributions have been made in 
the development of computational methods for solving Che equations of 
compressible vist'.ous flow. Chief among these has been the development 
of noniterative, bloch-trldiagonal Implicit Bcthods. These methods, 
which are not subject to restrictive stability conditions, are much 
more efficient than the earlier explicit methods. The newer methods are, 
however, BUch more complicated than the earlier ones and frequently still 
require long computation times. The goal of the present research is to 
develop a method for solving the compressible form of the Navier-Stokes 
equations at high Reynolds nueber that is unconditionally stable, com- 
putationally more efficient than existing methods, and simple and 
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•traifhtforvard to prograa. Th« Mthod to bo doocrlbod lo tht lapllclt 
onologue of tht oxpllclt flAlto-dlfforont aothod tht outhor prtotntcd 
In 1969 (ott (4J). Tht now aothod uott tht 1969 atthod ot Itt first 
stogt. Tht otcond ttogt rtaovto tht rostrlctlvt otoblllty condition of 
the 1969 atthod by rtcoitlng tht difftrtnce tquotlont in on iaplicit 
fom. Tht resulting aotrix tquatioos to bt solved srt either upper or 
lower block-bid isgonal equations and are aolvcd aorc easily than the 
block-trldiagonal aatrix equations of existing aethods. The acthod is 
second-order accurate in space and tiae and is presented in conservation 
form in two dimensions , Its extension to three dimensions is 
straightforward • 


II. The Navier-Stoxes Equations 

In two dimensions and by neglecting body force and heat sources, the 
unsteady compressible form of the Navier-Stokes equations may be written 
in conservation form as 


where 
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In terms of density p; x and y velocity components u and v; viscosity 
casifficients X and y; total energy per unit volume e; coefficient of 
heat conductivity k; and temperature T. Finally, the pressure p is 
related to the specific internal energy e and p by an equation of 
state» p(c, p), where e • e/p - (u^ + v*//2. 

The Navier-Stokes equations adequately describe aerodynamic flow at 
standard temperatures and pressures. If we could efficiently solve 
these equations there would be no need for experimental tests when 
designing flight vehicles or other aerodynamic devices. As John 
Von Neissann said in 19A6 [5], "Indeed to a great extent, experimentation 
in fluid dynamics is carried out under conditions where the underlying 
physical principles are not in doubt, where the quantities to be observed 
are completely determined by known equations. The purpose of experiment 
is not to verify a theory but to replace a computation from an unquestioned 
theory by direct measurements. Thus wind tunnels, for example, are used 
at present, at least In part, as computing devices to integrate the 
partial differential equations of fluid dynamics." 
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UoforfcuMttly th« •olutlon of tltoot oquotiono for flows of high loynolds 
BVAbors with strong vlseous-invlscld intsrsctlons hss dsfiod asthsMticsl 
•nslysls. Of tht two hoy fssturss of such flow - sspsrstlon and 
turbulsncs — vt have boon sblt to Bsks subscsntisl progrsss during the 
lest decade in the calculation of laoinar separation using maserical 
■ethoda. The calculation of turbulence largely raaains an unsolved 
probleo. Although the Havier-Stokes equations adequately describe such 
flows, computer speed and aemory llMtatlons make it impoaislble for the 
computational mesh to be fine enough in all spatial directions to resolve 
all significant eddy length scales of a hlgh-Reynolds-number turbulent 
flow. As Bradshaw said in 1972 16], "In turbulence studies we are 
fortunate in having a complete set of equations, the Navier-Stokes 
equations, whose ability to describe the motion of air at temperatures 
and pressures near atmospheric is not seriously in doubt (it is easy to 
show that the smallest significant eddies are many times larger than a 
molecular mean free path). We are unfortunate because numerical solu- 
tion of the full tlme^dependent equations for turbulent flow is not 
practical with present computers." 

In an approach that circiaovents the turbulence problem, the Reynolds or 
"time-averaged" Mavier-Stokes equations are solved. Thus, instead of 
seeking a time and spacially resolved solution of a rapidly fluctuating 
turbulent flow, only the time-averaged or mean flow solution is sought. 
This solution is sufficient to determine the principal quantities of 
interest, such as lift, drag, and heat transfer. The time-averaged 
equations look very similar to the original Navicr-Stokes equations 
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ftxctpc that •(]■« new ttraa* eaUad Rtynolda atraaa and turbulant haat- 
tranflar tania, appaar. Thaaa naw taraa rapraaant tha additional alxlng 
cauaad by turbulanea and ara datanlnad by aodala. Tha aodala vary Iron 
alapla algabraic axpraaalona to aata of additional dlffarantial oquationi 
that naed to ba aolvad. Although mich prograts haa alraady baan nada in 
the underatanding and modallng of tha phyaica of turbulanea, nueh nor a 
la needed before we will have tha capability to nunerlcally pradlct 
turbulent flow aeparatlon with confidence. 


HI. Humeric.il Method Applied to a Model Equation 

Before diacusalng the ntaoerlcal aolutlon of the Navier-Stokea aquatlona 
It la worthwhile to consider the aolutlon of the following alnpler nodal 
equation 


H ' It * 

The flow variable U govemerf by this equation convecta with speed c 
and diffuses with kinematic viscosity v. The laipllcit analogue of the 
author’s 1969 method applied to solve this equation yields the follow- 
ing set of finite-difference aquations [7]: 
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vhcre A Is shostn so that A > MX (|c| 4*(2v/Ax) - (AxMt)» 0.0). 


The above procedure contains tvo steps. The first step predicts a nev 
aoluclon at tine t ■ (n41)At at each nesh point 1 frost the known 
solution at time t ■ nAt, using a one-sided difference to approximate 
the first derivative term and a centered difference for the eecond 
derivative. The second stage of the predictor step enables the locally 
calculated solution changes to travel and diffuse throughout the 

flow field and then calculates Implicitly the solution change 
to be used In the third and final stage to determine the predicted 
solution The second step, or the corrector step* of the proce- 

dure Is similar except that It uses opposite one-sided differences to 
approximate first derivatives. 


The second stage of each step represents an Implicit approximation to 
the following equation 
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This «qu«cion daaerlbtt tha apraadlng of tha aolutlon ehanga Ac(au/»c), 
a tara of ordtr At» vlth apatd -X In tha pradictor and 4-X In tha 
corractor, Tha nat affact* if X ia of tha ordar of unity, la tha 
addition to tha aquationa of aotion of a tan of aacond ordar (tha 
differanca of two firat-ordar tana). Tha apraadlng aquation ia alao 
ralated to that obtainad by diffarautiating tha nodal aquation by t. 

The phlloiophy behind the procedure ia aa follova. Firat tha rata of 
Bolution change ia calculated locally at each naah point, uaing an 
explicit approxiaation to the governing phyaical aquationa. Thla local 
rate of change la only valid for a abort tine, equal approxinately to 
the tine required for a flow diaturbanca to travel froei one naah point 
to ita neighbor. Explicit procaduraa arc rcatrictcd, uaually for ata- 
bllity reaaons, to cioe atepi At laaa than or equal to thla charac- 
teristic disturbance transit tine. Second, this tine-step restriction 
is rraoved in the second stage by allowing the locally daterained rates 
of solution change to convect and diffuse globally throughout the flow 
field, governed by an equation related to the phyalca of the flow. This 
latter equation is solved Inplicitly to dateminc the rate of solution 
change at each point that approxinatea the actual rata during the 
entire interval At. 

The nethod is, according to linear theory, unconditionally atable 
(unbounded At), requires the solution of bidiagonal aquations only, 
and ia second-order accurate vndcr the constraint that vAt/Ax* 
raaalns bounded as At and Ax approach aero (i.c., X ranains of 
the order of unity). 
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Mot* that if th« quantity k ia larof tba aacond atagaa raduea ta 

♦1^ • A»J 


«bJ*‘ - 41^* 

or no iapllclt proeadura at alli and tha aathod ii idantical to tba 
1969 narhod. Such a cholca for k raaulta if tha choaan tlaa atap it 
alraady aatiafiaa tha acability condition of tha 196P Mthod, 


Bacauae of thic faaturc, the nathod haa an advantage in mnerical affi- 
Clancy ovar axlatlng Inplicit aathoda. Not only ara the ntaarical pro- 
cedures alnpler — bidlagonal versus tridiagonal — but in flou regions 
for which it aatisfles the a?»nve liability condition! the nethod 
reduc'is from an iaplicit to a ainplcr explicit onCi 


IV, Nuawrical Method Applied to the Navier-Stokea Equations 
Applying the nethod to aolve the Navier-Stokea aquations we obtain the 
following iaplicit predictor-corrector sat of finite-difference 
aquations. 
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where A^/Ax , A^/Ax, Aj^Ay and A__/Ay /iira dlffcrance operators daflnid by 
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‘-‘l.3 - ‘l.3 ~ ^1.3-> 

Ay Ay 


As for the model equation, the first derivative terms are one-sided 
differenced (as shown above) and the second derivative terms arc centrally 
differenced. The matrices |a| and |b| arc matrices with positive eigen- 
values and are related to the Jacobians of F and G. Let S , S , and 

X y 

their inverses denote the matrices that diagonalise A and B with 
M > 1 • k ■ 0 (viscous terms neglected). If the gas equation of state 
is perfect, p ■ (y - l)pc, A • ® " ®y^^B®y’ 
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ind whir® c ■ VYf/P li the epeed of aound, a • (1/2) (u* + v*) and 
0 « Y - 1* 

The netriceB S and S are each given above aa the product of two 

" y 

matrices, For each* the right matrix repreaenta a transformation from 
conservative to nonconservative variables, for example, from (6p, 6pu, 
6pv, «e) to (dp, 6u, dv, 6p). The left matrix transforms from noncon- 
servative to characteristic form (dp - dp/c*, pcdu -f dp, dv, -pdu dp) 
and (dp - dp/c®, du, pcdv + dp, -pcdv “f dp) for the and utrices, 
respectively. The inverses s”^ and s“^ are simple to derive. The 
matrices |a| and |lij arc defined by 

1*1 • ““ 1*1 • 
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V > sax (Mi X 2 m, k) 


ViscouB effects are included through the use of the viscous coefficient 
V. For some test probleos this coefficient had to be increased during 
the initial part of the calculation when large transients in the solu- 
tion occurred. 
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For rafioM of tha f,low in which At aatiafiaa tha following axplieit 
atability eonditiona, 

gt < 1 ^ - — 

* (|w| + c)/Aic 't (2v/pAx*) 
and 

^ I I 

" ^ (|v| + c)/Ay (2v/pAy*) 

all and all X, vanlah and tha oat of diffarance aquations raducas 
A s » 

to tha 1969 explicit aquations. For other regions In which neither 
relation Is satisfied, the resulting difference equations are either 
upper or lower block bldiagonal equations with fairly straightforward 
solutions. 

V. Numerical Results 

The method was applied to solve for the Interaction of a ahock wave 
litcident upon a boundary layer. The flow field is sketched in Fig. 1. 

As shown In Fig. 2, the initial coniitlor. was that of uniform flow, and 
the condition at the top mesh boundary W«s such that a shock wave of given 
strength would be generated and impinge upon a flat plate at a given 
point. The conditions at the upstream boundary were held fixed at their 
initial supersonic free-stream values; the downstream boundary conditions 
were obtained by sero-order interpolation; the lower boundary conditions 
were obtained by reflection. 

The mesh contained 32 H 32 points, with 16 spanning the boundary layer. 

The time step was chosen so that the free stream moved approximately IX 
during each time step. With this choice tha time atep satisfied the 
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•xpliclt Btability condieioiiB •vcrywhBr* 4Mc«pfc In the fint M«h Bpannlni 
ttkc boundary layer. 

In Fig. 3| the coaputad reaulta arc coaparad with axpcrlacnt (8) and 
vlth boundary-layer theory in the abacnce of a ahock wave (9}. Ihc reaulta 
are for Mach 2 laalnar flw at a Rcynolda niaibcr of 2.95 * 10*. The 
calculation (1) uaed Suthcrland'a foraula to calculate aolccular viccoalty; 
(2) vaa run for 256 tine ateps, at which tine the aeah vaa rezoned to 
cover just the interaction region; and (3) was run for an additional 
256 tine steps. It required about 1.5 nin of conputer tine on a CI)C 7600 
conputer. The results for skin friction and surface pressure compare 
favorably with those of theory and experiment. 

In Fig. 4 the calculated velocity profiles ahead, within, and aft of the 
separation region are compared with the conputad results, using the 1969 
nethod alone. The two sets of results agree closely; however, the com- 
puter time required by the newer method was more than an order of magni- 
tude less than that required by the 1969 method. 

The computation tines for a aeries of laminar and turbulent boundary- layer 
interactions with shock waves are given in Table 1. For each problem the 
flow was computed to the aame physical time, which for the new method 
required 256 time steps. For the turbulent flow cases, a aimpie algebraic 
eddy viscosity model {10] was used to account for the effects of turbulence. 
For each case, the table shows the Reynolds number, the ratio of the 
time step uaed to the maximum allowed by explicit stability conditions 
(Courant-Friedrich-Lewy number, or CFL), the computer time required per 
time step per grid point on a CDC 7600 computer, and the total computer 
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tlM raqulrtd. The tabulated raaulta ahow that tha naw aathod la ona to 
thraa ordara of aagnituda aora afflclant than tha 1969 aathod. 

For tha taat Ci.aaa eonaldaredi the nav aathod la aora than tvica aa faat 
par tiaa atap par grid point an tha bloek-tridiagonal implicit aathoda 
in uaa today. Fart of tha raaaon for thin ia tha aaih-poirtt apacing and 
tiaa atap choaan; of tha total nuabar of aaah pointa aora than half 
raquirad only tha uaa of the 1969 explicit aathod. At theaa aaah pointa 
tha choaan time atap already aatiafiad the local explicit atabllity con- 
ditioni therafore, the implicit procedure, the aecond atage, waa omitted, 
Tha implicit procedures were raquirad only in the fine maah apannlng the 
boundary layer, where explicit stability conditions would have imposed a 
severe time-step restriction. It is estimated that if the Implicit pro- 
cedures were used at each grid point, the time step per grid point for 
a two-diatensional calculation would be 2.45 x lO"** sec for laminar flow 
and 2.75 x lO”** sec for turbulent flow. The difference between these last 
two values represents the additional computation needed to evaluate the 
tuT ulence model relations. 

VI. Conclusion 

A new ntaserical method has been devised for solving the equations of 
compressible viscous flow. The method represents the Implicit analogue 
of the explicit method presented by the author in 1969. It is uncondi- 
tionally stable, second-order accurate, and, for many applications, is 
more efficient than other methods in use today. Because the new method 
uses the 1969 method as its first stage, many existing computer programs 
in which the 1969 method is used can be updated by adding the deacribed 
implicit aacond atage. 
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Tablt 1 Coaputation tiat 


Cat* 

Method 

CFL 

CDC 7600 tlaa 
atap grid point 

Total 

tiaa 

Laalnar 

1969 

0.9 

1.25 K 10"^ aac 

12 ain 

R ■ 3 K 10* 

Maw 

20 

1.55 * 10“* aac 

41 aac 

Turbulen'c 

1969 

.9 

1.55 >« 10’^ aac 

2 hr* 

R • 3 * 10‘ 

Maw 

160 

1.85 X 10"" aac 

48 aac 

Turbulent 

1969 

.9 

1.55 X 10"" aac 

15 hr* 

R - 3 X lO’ 

New 

1200 

1.85 X 10"" aac 

4B tec 


*E8tiauit«d, 
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ri«urt Caption! 


Fig. 1 Sketch of ohock. boundary-flayer interaction 

Fig. 2 Initial flow field for ahock, boundary- layer interaction 

Fig. 3 Compariion of raaulta. a) Surface praaaure. b) Skin friction 

Fig. 4 Conpariion of velocity profilaa 
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